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ABSTRACT 

Understanding the chemical evolution in star-forming cores is a necessary pre-condition to correctly 
assess physical conditions when using molecular emission. We follow the evolution of chemistry and 
molecular line profiles through the entire star formation process, including a self-consistent treatment 
of dynamics, dust continuum radiative transfer, gas energetics, chemistry, molecular excitation, and 
line radiative transfer. In particular, the chemical code follows a gas parcel as it falls toward the center, 
passing through regimes of density, dust temperature, and gas temperature that are changing both 
because of the motion of the parcel and the evolving luminosity of the central source. We combine 
a sequence of Bonnor-Ebert spheres and the inside-out collapse model to describe dynamics from 
the pre-protostellar stage to later stages. The overall structures of abundance profiles show complex 
behavior that can be understood as interactions between freeze-out and evaporation of molecules. 
We find that the presence or absence of gas-phase CO has a tremendous effect on the less abundant 
species. In addition, the ambient radiation field and the grain properties have important effects 
on the chemical evolution, and the variations in abundance have strong effects on the predicted 
emission line profiles. Multi-transition and multi-position observations are necessary to constrain the 
parameters and interpret observations correctly in terms of physical conditions. Good spatial and 
spectral resolution is also important in distinguishing evolutionary stages. 
Subject headings: 



1. INTRODUCTION 

In star-forming regions, the chemistry is not in steady 
state while cores are being formed by contraction of parts 
of molecular clouds, or later, while the cores are collaps- 
ing. Since star-forming cores can be studied through 
emission and absorption by gaseous molecules, it is nec- 
essary to understand the chemical evolution in order to 
correctly assess the physical conditions. Dust continuum 
emission is capable of providing information on the core's 
physical structure and evolutionary state (Myers & Ladd 
1993, Andre, Ward-Thompson, & Barsony 1993), but ve- 
locity structure in the core, the other important physical 
parameter to constrain theoretical models in star for- 
mation, can only be obtained from molecular spectra 
(Mizuno et al. 1990; Zhou et al. 1993; Choi et al. 1995; 
Gregersen et al. 1997; Mardones et al. 1997; Gregersen 
et al. 1998; Di Francesco et al. 2001; Belloche et al. 
2002). However, the observations of molecular spectra 
can be misinterpreted if the molecular abundance distri- 
butions are not well understood. For example, Rawings 
& Yates (2001) caution that the blue asymmetry, which 
has been used to detect infall in star-forming cores, could 
disappear by the depletion of molecules in spite of infall 
velocity structure. Therefore, we need to understand the 
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chemical evolution during star formation to isolate good 
molecular tracers for each evolutionary stage and differ- 
ent regions along the line of sight. 

However, chemical processes in star-forming regions, 
especially reactions on grain surfaces and the gas-dust 
interaction, are poorly understood. Several theoretical 
studies (Bergin & Langer 1997; Aikawa et al. 2001; 
Caselli et al. 2002; Li et al. 2002; Shematovich et al. 
2003; Aikawa et al. 2003) examined the chemistry in pre- 
protostellar cores (PPCs). These are considered as the 
potential sites of future star formation because they are 
believed to be gravitationally bound, but have no central 
hydrostatic protostar (Ward-Thompson, Scott, & Andre 
1994; Ward-Thompson, Motte, & Andre 1999). The only 
heating source of the PPCs is the interstellar radiation 
field (ISRF), so the dust temperature in the inner denser 
regions is as low ~7 K (Evans et al. 2001; Zucconi, 
Walmsley, & Galli 2001). Those theoretical studies men- 
tioned above modeled the gas-phase chemical evolution, 
including gas-dust interactions, and followed the chang- 
ing depletions of some molecules and ionization fraction 
with time/density. Each assumed that the gas tempera- 
ture {Tk) is equal to the dust temperature {Td) and con- 
stant through the cloud. The consistent result from those 
studies is that sulfur- or carbon-bearing molecules are 
easily depleted from the gas in dense inner regions dur- 
ing early stages, but nitrogen-bearing molecules become 
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abundant at later times without significant depletion. 
Based on the conclusion of those studies, we can consider 
nitrogen-bearing molecules such as N2H+ as good tracers 
of dense gas in evolved PPCs. However, N2H"'" may also 
deplete (Bergin et al. 2002) in PPCs due to the long dy- 
namical timescale compared to depletion timescale. In 
addition, Lee et al. (2003) showed that the chemical 
evolution may depend on the dynamical scenario (e.g. 
free-fall or ambipolar diffusion) or different properties of 
dust grains as well as different density structures. 

In contrast to the situation in PPCs, there are few 
studies of chemistry in cores after collapse begins. This 
is because the chemistry becomes more complicated due 
to additional heating mechanisms such as accretion lumi- 
nosity or the luminosity from protostellar objects. As a 
result of the additional heating sources, the dust temper- 
ature increases, and icy materials will be evaporated from 
grain mantles, changing the available molecular probes. 
Rawlings & Yates (2001) tried to calculate the evolution 
of chemistry and line profiles in Class sources using 
a self-consistent chemical and dynamical model. They 
adopted the multi-point chemical model of Rawlings et 
al. (1992) in order to constrain the observed line pro- 
files. However, they used only one dust temperature pro- 
file, obtained from dust continuum observations of B335 
(Zhou et al. 1990), assuming that the temperature profile 
docs not vary with time. In addition, they did not include 
gas energetics. Ceccarelli, HoUenbach, & Tielens (1996) 
made the first attempt to calculate self-consistently line 
spectra from the infalling gas at wavelengths between 
20 /zm to 200 /zm with a model that includes dynamics, 
chemistry, heating and cooling, and line radiative trans- 
fer. In this work, they separated the dynamical and ther- 
mal evolution of the model system. For dynamics, they 
adopted the simple solution of spherical isothermal col- 
lapse (Shu 1977). The dust temperature was taken from 
the solution of Adams & Shu (1985) without calculat- 
ing dust continuum radiative transfer. In addition, they 
only considered 01 fine-structure lines and CO and H2O 
pure rotational lines as cooling mechanism, even though 
they included NIR pump heating and dynamical com- 
pressional heating, which are important heating sources 
at the inner hotter region, as well as gas-grain collisional 
heating. According to their results, the difFerence be- 
tween dust and gas temperatures is ~25% for the inner 
denser region, but the difference is larger in the outer 
layers. Other heating processes which could significantly 
contribute to the gas energetics, such as photoelectric 
heating and cosmic-ray heating, were not considered by 
Ceccarelli et al. For chemistry, they considered 182 pure 
gas-phase reactions including only 44 species, most be- 
ing oxygen- and/or carbon-bearing molecules. In their 
work, water vapor, at the inner region where Td > 100 
K, evaporated into the gas phase so that the cooling by 
H2O rotational lines was enhanced significantly in the hot 
region. The most inconsistent assumption in this work 
with the results of studies for PPCs was the constant ini- 
tial abundance profile of molecules. They assumed the 
timescale of evolution before collapse was much shorter 
than the chemical timescale, so the chemical composition 
remained unchanged during the first phase of the con- 
densation of the matter into a singular sphere. However, 
we know from observational studies and the theoretical 
models described above that the abundance profiles in 



PPCs are not constant. 

More recently, Rodgcrs & Charley (2003) combined the 
dynamical evolution with the chemical evolution focus- 
ing on the hot core phase, where ices evaporate from 
grain surfaces, and high-temperature chemistry is im- 
portant. They showed that the dynamical timescale in 
Shu's inside-out collapse model was too short compared 
to the chemical timescale to develop the observed molec- 
ular abundances in massive hot cores. Doty et al. (2004) 
also constructed a self-consistent radiative transfer model 
including dust radiative transfer, gas energetics, and a 
time-dependent chemical network to simulate molecular 
line transitions and, finally, to compare with observations 
in IRAS 16293-2422. They adopted a static protostellar 
envelope and assumed initial abundances based on ob- 
servations of massive star-forming cores. They did not 
include the freeze-out process in the chemistry and con- 
sidered that only CO, H2CO, and CH3OH were desorbed 
in the outer envelope that has temperatures below 100 
K. 

In addition to these theoretical studies, Lee et al. 
(2003) and J0rgensen et al. (2004b) have tested some 
simple empirical abundance profiles, such as a step func- 
tion and a drop function (see Figure 9 of J0rgensen et al. 
2004b) to fit observed line profiles of some molecules in 
the PPC stage and in the Class stage, respectively. 

In this paper, we calculate the evolution of chemi- 
cal abundances and line profiles of several important 
molecules, which have been used to study star-forming 
cores at submillimeter wavelengths, in a more self- 
consistent fashion than seen in earlier work. We include 
a specific dynamical model, the dust continuum radia- 
tive transfer, the gas energetics, a chemical model, and 
the line radiative transfer. Our models include evolution 
in the PPC stage and the later evolutionary phase after 
collapse begins. 

In §2, we describe all models used for this study, and 
§3 and §4 show results of our modeling. We discuss in §5 
how various molecular line tracers depend on evolution- 
ary stage, cloud environment, etc. Finally, we summarize 
our conclusions in §6. 

2. MODELS 

In this section, we describe our modeling process from 
dynamics to line profile simulation through continuum 
radiative transfer, gas energetics, chemical evolution, and 
line radiative transfer. The flow chart in Figure 1 shows 
the sequence of the modeling. First, density and veloc- 
ity fields [n(r) and u{r)] as a fimction of time are pro- 
vided by star formation theory. Most dynamical theories 
do not provide the temperature distribution [T(r)], but 
theories do usually provide a way to calculate the lumi- 
nosity of the central object versus time [i(i)]- There- 
fore, we can calculate the dust temperature distribution 
[TD(r)] in the envelope at any point in the evolution by 
calculating the energy transfer through the dust. In this 
dust continuum radiative transfer calculation, we have 
to assume the opacity of dust grains («;,y) and the gas 
to dust mass ratio as well as calculating the luminosity 
of the central object. Once we calculate the dust tem- 
perature distribution, the gas temperature distribution 
[Tff(r)] can be computed from the gas energetics. Deep 
in the cloud, where densities are high, collisions between 
gas and dust maintain Tk = Td, but, at lower densities. 
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Tfc must be calculated by balancing heating and cooling 
rates. In order to simulate molecular line profiles, we 
need to know the distributions of molecular abundances, 
so we adapt the chemical model of Bergin et al. (1995) 
to compute the chemical abundances as a function of ra- 
dius [-'^(r)] at any evolutionary time. We calculate the 
abundance changes following every parcel of a model core 
through the whole evolutionary time in a Lagrangian co- 
ordinate system because the abundances established in 
the outer layers are carried inward with the matter. Af- 
ter the calculation, all parcels are transferred to Eulerian 
coordinates so as to give the distributions of molecular 
abundance at any given time. The next step is to calcu- 
late level populations of molecules with the Monte Carlo 
(MC) method to solve the coupled excitation and the 
radiative transfer self-consistently. The final step is to 
simulate line profiles at each time step by calculating the 
intensity emitted by the model core at the distance of d 
over the spectral line velocities at any impact parameter, 
convolving to match the spectral {6v) and spatial reso- 
lution {9pwhm), and the main beam efficiency (rjmb) of 
the observations. Each step of the whole process of this 
modeling is described in more detail in following subsec- 
tions. 

2.1. Dynamical Model 

We use the inside-out collapse model (Shu 1977) as a 
test dynamical collapse model. However, we need to con- 
sider the initial conditions in the chemistry carefully be- 
cause previous studies have shown that many molecules 
are depleted in dense cores (PPCs), whose central density 
is high enough to make a big contrast to their boundary 
density, but which do not contain a source of central lu- 
minosity. 

The submillimeter continuum emission from those 
PPCs indicates that the density profile does not follow 
a power law all the way to the center. Instead, it shows 
roughly a flat density distribution within some radius. 
Such density profiles are well described by Bonnor-Ebert 
spheres. Evans et al. (2001) showed that Bonnor-Ebert 
spheres fit the observed data well through detailed analy- 
sis of the submillimeter emission from dust toward three 
PPCs (sec also Alves, Lada, & Lada 2001). The cores, 
whose only heating source is the interstellar radiation 
field (ISRF), are very cold in the interiors (To ~7 K). 
As a result, the densities needed to match the observed 
emission arc higher than what had been thought previ- 
ously. This is important for studies of molecular evolu- 
tion because higher densities enhance the freeze-out rate. 

For higher central densities, the region of the flat den- 
sity distribution in Bonnor-Ebert spheres is smaller, and 
the outer regions are similar to a power law (n(r) oc r"^) 
with p ~ 2. As the central condensation is growing, 
the spheres approach the singular state, which is used 
for the initial condition of inside-out collapse. There- 
fore, we combine a sequence of Bonnor-Ebert spheres 
and inside-out collapse in order to set up a dynamical 
model to evolve from the pre-protostellar stage to later 
stages. The outer radius for both Bonnor-Ebert spheres 
and inside-out collapse used in this work is 0.15 pc. We 
assume a constant density structure of 10^ cm~^ as the 
very first condition, and the model core is embedded 
in surroundings with some amount of visual extinction. 
Figure 2 shows the density profiles of 7 Bonnor-Ebert 



spheres and the density profiles of inside-out collapse in 
several time steps. The (7) Bonnor-Ebert evolutionary 
steps assume an evolution in central density from 10** to 
10^ cm~^. The Bonnor-Ebert spheres do not evolve dy- 
namically, so there is no good constraint on timescales 
in the PPC stage. Therefore, we do not know the exact 
timcscalc spent at each central density, but can expect 
that the timcscalc decreases as the central density in- 
creases. Based on the free-fall timcscalc (~ l/^/n) and 
the statistical timescale of the whole PPC stage (see e.g. 
Lee & Myers 1999), we simply assume the timescales of 
2.5 X 10^, 1.25 X 10^ 6 x 10'', 3 x 10'*, 2 x 10^, 1 x 10'', and 
5 X 10"^ for Bonnor-Ebert spheres of ric ~ 10'', 3 x 10'', 
10^ 3x 10^ 10*^, 3x 10^ and 10^ cm^^ respectively. The 
timescale for the very initial constant density (n = 10^ 
cm~^) is 5 X 10^ years. Therefore, the total timescale for 
the PPC stage is a million years. We, however, expect 
to give some constraints to the timescales by comparing 
the results of this work with actual observations. For 
instance, we know that the most important effect of the 
PPC stage is the freeze-out of molecules, so the longer 
this stage lasts, the more the lower density in the outer 
layers will be affected. 

We obtain density profiles from the Bonnor-Ebert 
sphere model and inside-out collapse model assuming 
that the model core is isothermal with the kinetic tem- 
perature of 10 K. Evans et al. (2001) showed that differ- 
ences between the isothermal Bonnor-Ebert spheres and 
non-isothermal Bonnor-Ebert spheres are small. The to- 
tal core mass is 3.6 Mq. The properties of the model 
core are summarized in Table 1. 

In Bonnor-Ebert spheres, all gas parcels are in steady 
state, but the density at a given position increases as time 
goes by without any noticeable flow of gas. In collapse, 
all gas parcels inside the infall radius move inward with- 
out any interaction with adjacent gas parcels, but gas 
parcels outside the infall radius stay at the original radii 
without any change of density distribution with time. In 
order to study the chemical evolution in a model core, 
we should follow every gas parcel that experiences dif- 
ferent physical environments with time (best done with 
a Lagrangian coordinate system) and calculate chemi- 
cal reactions adopting the results of chemical evolution 
from one time step earlier as the initial condition at a 
given time step. Finally, we transfer all gas parcels again 
to Eulerian coordinates for each given time step to de- 
termine the abimdance profiles of molecules with time. 
We consider 190 time steps in the inside-out collapse to 
make total 198 time steps for the whole evolution. The 
time step varies with time. The velocity change slows 
down with time, so a smaller time step has been used for 
earlier times; 500 years for the first 10,000 years, 1000 
years from 10,000 to 100,000 years, and 5000 years from 
100,000 to 500,000 years. 

We use the continuity equation in Lagrangian coor- 
dinates to realize this idea in practice. The continuity 
equation in spherical coordinates is 
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where 1? is a velocity vector, and u is the radial compo- 
nent of velocity. The change of density of a parcel in a 
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small time interval, dt is given by 

,2u du, , 

dn = -n\ h -7^]dt. (2) 

r or 

Therefore, the density of a parcel in the next time step 
is given by 

n{t2) = n{ti)+dn = n{ti)-n{ti)[—+—]dt, dt = t2-ti 

(3) 

where ti and t2 represent a given time step and the next 

time step, respectively, and u and ^ is the mean ve- 
locity and mean velocity gradient between ti and t2, 
respectively. We use mean velocity and mean velocity 
gradient because the velocity in collapse changes contin- 
uously with time, and the time interval for this work is 
not infinitesimal. The other reason that we use the mean 
values is in the propagation of the infall radius outward 
at the speed of sound. The region between two infall 
radii at ti and t2 has zero speed at ti, so we cannot use 
u(ti) and (ii) to calculate the change of density of a 
parcel. We assume c?n = if dn has a negative value, dn 
can have a negative value between infall radii at ti and 
t2 because u{ti) = at those radii so that tl is a half of 
u{t2), and in turn, its absolute value can be smaller than 

that of which is always positive. 

Once we calculate n{t2), we can find the radius, tem- 
perature, and velocity of the given parcel from the solu- 
tions of inside-out collapse at t2- We follow all 512 gas 
parcels until they fall onto the central protostar. Wc as- 
sume that a gas parcel falls onto the central protostar if 
the radius of the parcel is less than the inner radius of 
the model. 70 AU. The falling time depends on the ini- 
tial position; that is, the smaller the radii gas parcels are 
at initially, the faster they disappear from the envelope. 
Figure 3 shows the evolution of density, temperature, 
radius, and Ay of several parcels with time. Ay is cal- 
culated from outside to a given parcel. Wc consider the 
initiation of collapse as i = 0, so the evolution is divided 
into t < and t > for the Bonnor-Ebert spheres and 
for the inside-out collapse, respectively. 

2.2. Dust Radiative Transfer 

In order to calculate dust temperature profiles (TD(r)), 
we use a 1-d dust continuum radiative transfer code, 
DUSTY, created by Ivezic et al. (1999). PPCs are 
heated only by the ISRF (Evans et al. 2001). For t > 0, 
we add a central luminous object. Even for t > 0, the 
ISRF significantly affects both the total observed lumi- 
nosity and the shape of the observed submillimeter in- 
tensity profile (Shirley et al. 2002, Young et al. 2003). 
Evans et al. (2001) found that observed data were better 
fitted with the ISRF reduced for wavelengths from the ul- 
traviolet to far-infrared. We assume that the model core 
is exposed to the Black-Draine ISRF (Evans et al. 2001) 
attenuated by some amount of external visual extinction 
(Ay^o) in the dust radiative transfer code because we 
test several external visual extinctions between 0.5 and 
3 mag in chemical evolution. Our standard model for 
the chemical evolution adopts Ay^e = 0.5 mag. For the 
opacity, wc use 0H5 dust, which is found in column (5) 
of Ossenkopf & Henning (1994). This opacity, for grains 
that have coagulated and accreted thin ice mantles, was 



found to match the multi-wavelength observations of low 
mass cores in all stages modeled here (Evans ct al. 2001, 
Shirley et al. 2002, Young et al. 2003). Figure 4 (left 
panel) shows T^ir) before collapse begins {t < 0). The 
dust temperature goes down to about 6 K at the core cen- 
ter because the model core is heated only by the ISRF 
at < < 0. 

For the evolution of the central source of luminosity 
(for t > 0), we adopt the model of Young & Evans (2004). 
After collapse begins, the central luminosity source in- 
cludes accretion from the envelope onto a central proto- 
star and disk, the disk luminosity, and the internal stel- 
lar luminosity due to contraction and deuterium burn- 
ing. Some of these sources turn on at different times. 
At early times, the accretion luminosity dominates all 
other components. Since there is no disk at very early 
times, accreting material falls directly onto the central 
object. The first hydrostatic core (FHSC) (Larson 1969, 
Boss 1993, Masunaga et al. 1998), which has a large 
radius (5 AU), is included until 20,000 years. The tran- 
sition between the large core to the smaller core occurs 
between 20,000 and 52,000 years. As the disk begins to 
form, most accreting material is deposited onto the disk, 
so the accretion luminosity becomes the dominant lumi- 
nosity source until the photospheric luminosity becomes 
important at t = 100, 000. 

Figure 4 (right panel) shows the dust temperature pro- 
files calculated from the dust radiative transfer code for 
the model core after dynamical collapse begins. In gen- 
eral, temperature increases with time at all radii as the 
central luminosity increases, and the envelope mass de- 
creases, but there are two jumps between 2x10** and 
5x10"^ and between 10^ and 1.6 x 10'"' years. The first 
jump is caused by the transition from the FHSC to the 
smaller core, and the second jump is due to the turn-on 
of the photospheric luminosity from gravitational con- 
traction and deuterium burning in the central star at 
t ^ 1 X 10'"' years. Around 100000 years, the transition 
between Class and Class I occurs, when the bolometric 
temperature is about 70 K. 

2.3. Gas Energetics 

The gas kinetic temperature (Tr-) is determined by the 
balance between heating and cooling. The gas heating 
at large densities is dominated by gas-grain collisions, 
but grain photoelectric (PE) heating is significant at the 
outer part of a core. Cosmic-ray (CR) heating is also 
important in a transition region where PE heating has 
decreased and gas-grain collisions are not yet dominant. 
Therefore, we use a new code updated by S. Doty to 
add PE and CR heating to the code of Doty & Neufeld 
(1997), which considered only gas-grain coUisional heat- 
ing and adopted the cooling functions of Neufeld, Lepp, 
& Melnick (1995) to determine the gas cooling rates due 
to emission from CO, O2, O I, H2, II2O as well as from 
other diatomic and polyatomic species. They computed 
the steady-state abundances of the atoms and molecules 
in pure gas chemistry including the full UMIST chemical 
network (Millar ct al. 1991). 

However, molecules such as CO (the dominant coolant 
in cold dark clouds) are frozen onto dust grain surfaces 
at very dense and cold inner region of PPCs. In addi- 
tion, once a central star forms dust grains can be heated 
above 100 K, where the evaporation of water from icy 
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grain mantles occurs. Water is a very efficient coolant 
in hot regions. Therefore, one should, in principle, in- 
clude the effects of freeze-out and desorption to compute 
the correct abundances of coolants in the dense inner re- 
gion, and iterate the calculation of heating and cooling 
to find the converged gas temperature for different densi- 
ties and dust temperatures. However, Goldsmith (2001) 
and Doty & Neufeld (1997) found the effects of depletion 
of CO and evaporation of water were relatively small at 
high densities because of efficient gas-dust coupling. We 
tested the effect of CO depletion by reducing X(CO) by 2 
orders of magnitude in the gas energetics code. The cal- 
culated gas temperatures are the same as the dust tem- 
peratures if the density is higher than about 10^ cm^^, 
where the freeze-out of CO becomes significant. In addi- 
tion, the inner radius in our models is not small enough 
to reach layers where the dust temperature is >100 K. 
Therefore, we do not include the effects of freeze-out and 
desorption in the gas energetics. Other possible inputs 
into the gas energetics are H2 formation, FUV pumping 
of H2, and H2 dissociation heating. In our calculations 
wc have assumed a weak attenuated UV field leading H2 
to be fully shielded well before the layers examined by 
our models. Therefore the heating associated with H2 
formation, dissociation, and FUV pumping will be small 
compared to photo-electric heating and can be ignored 
(sec Sternberg & Dalgarno 1989; Burton, Hollenbach, & 
Tielens 1990). 

The gas energetics code can examine solutions with 
an attenuated UV field. The CO lines, which are easily 
thermalized, around PPCs show much weaker emission 
lines than those expected if the cores are heated by the 
unshielded local ISRF (K. Young et al. 2004). In order 
to calculate Gq, which is the strength of the UV field 
relative to the local ISRF at the outer cloud boundary, 
we use the equation. Go = exp{—1.8 x Ay.e)- Figure 5 
shows gas temperature profiles computed with this gas 
energetics code. The gas temperature gets high in the 
outer regions because of photoelectric heating. However, 
gas-grain collisions are dominant in dense regions, so the 
gas temperature is well coupled with the dust tempera- 
ture at small radii. If the ISRF suffers less attenuation, 
then it heats gas better through the core to keep the gas 
temperature higher, especially, in less dense regions. 

Details of the gas energetics code used in this study 
are given in the appendix of K. Young et al. (2004) and 
Doty & Neufeld (1997). 

2.4. Chemical Model 

We use the time-dependent chemical model of Bergin, 
Langer, & Goldsmith (1995), updated in Bergin & 
Langer (1997), which includes depletion and desorption 
of atoms and molecules onto and off grain surfaces, as 
well as gas-phase interactions. However, the model does 
not explicitly include chemical reactions on grain man- 
tles. To reduce the computation time we have used a 
reduced network including ^80 species and ^800 reac- 
tions. The reduced network was created by an exami- 
nation of the major formation and destruction paths for 
key species. This was tested in comparison to the larger 
network. 

As described in §2.1, chemistry is calculated following 
each parcel, which experiences different density and tem- 
perature environments as it falls toward the center. We 



consider 198 dynamical time steps from the PPC stage 
to later stages after collapse begins as mentioned in §2.1. 
In each time step, physical conditions are constant, so 
the chemistry is calculated pseudo-time-dependently be- 
tween two time steps by adopting the chemical results 
of the previous time step as the initial conditions in the 
later time step without changing physical conditions. 

This model for gas-phase chemistry includes the pho- 
todissociation and photoionization of molecules as well 
as the cosmic ray ionization. For the cosmic ray ioniza- 
tion rate we have assumed Ci?2 = 1-2 x 10~^^ s~^ and 
the abimdancc of He"*", a key destroyer of many species, 
is derived from a balance of cosmic ray ionization of He 
(Ci/e+ — 6.5 xlO"-^^ s~-^) and its major reactants (espe- 
cially CO). However, to reduce further the computational 
complexity we have not included CO photodissociation 
in the model and only allow its abundance to be altered 
via gas-phase and gas- grain interactions. To minimize 
the effects of this on the chemical calculation we have 
assumed there exists additional extinction (Av,e) beyond 
the edges of our model. CO is self shielding for Av be- 
tween 0.5 and 1 mag, so most of our results arc unaffected 
by neglecting CO photodissociation. In our models we 
have examined external extinctions between 0.5 - 3 mag. 
At its lowest value the abundances of simple radicals (e.g. 
CN, C2H) that form in the outer photo-dissociation lay- 
ers could be affected. 

For the gas-grain reactions, this model assume that 
grains are charged with one electron per grain, and pos- 
itive ions undergo dissociative recombination with the 
same branching ratio as in the gas phase when they col- 
lide with grains. This model also considers three im- 
portant desorption mechanisms: thermal evaporation, 
cosmic-ray-induced heating, and direct photo-desorption 
(see Bergin, Langer, & Goldsmith 1995 for details). At 
low temperatures as in PPCs, cosmic-ray-induced des- 
orption is dominant. Once collapse begins, however, Tjj 
increases, so thermal evaporation becomes important. 
Freeze-out of molecules is highly dependent on the dom- 
inant component of the grain mantle, that is, the surface 
binding energy of the species. We use the surface bind- 
ing energy of species onto a bare Si02 dust grain as a 
standard. Table 2 lists the binding energies of selected 
species. We also test binding energies onto H2O and CO 
dominant grain mantles. We assume a sticking coefficient 
of 1.0 for all species. In order to reduce the water abun- 
dance in the gas to match SWAS observations (Bergin et 
al. 2000), we adopt a binding energy of 5700 K (Fraser 
et al. 2001) for water onto grain surfaces. This binding 
energy is the same even on different dust grains, unlike 
the binding energies of other molecules. In addition, we 
consider that atomic oxygen is frozen onto grain surfaces 
in the form H2O on grains again to match SWAS obser- 
vations. For photodesorption, this model uses the mean 
interstellar UV field (Habing 1968), which is attenuated 
as a parcel moves inward based on the calculated Ay- 
We have used the Black-Draine field for the dust heating 
calculations and the Habing field for all calculations of 
photodissociation and photodesorption. This small, fac- 
tor of 1.7, inconsistency will not affect our results given 
uncertainties in photodissociation rates and photodes- 
orption yields. We assume that the core evolution occurs 
in a shielded region, with different external extinctions 
to represent various environments of star- forming cores. 
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The initial elemental abundances are listed in Table 3. 
Initially (t= —10^ years), all hydrogen is molecular (H2), 
and all carbon is ionic (C"*"). 

For more detailed explanation of this chemical model, 
refer to Bergin et al. (1995). 

2.5. Molecular Radiative Transfer Model 

Once all physical conditions [n(r), TKir), u{r), X{r)] 
are computed from other models we perform a radiative 
transfer calculation for the spectral lines. For the molec- 
ular radiative transfer calculation, we use the MC code 
that has been developed by Choi ct al. (1995). The 
MC code generates model photons at a random position, 
in a random direction, and at a random frequency with 
proper random number distributions. The MC code cal- 
culates the excitation by the model photons and uses sta- 
tistical equilibrium to adjust each level population until 
the criteria of convergence arc satisfied. We can model 
arbitrary distributions of systematic velocity, density, ki- 
netic temperature, microturbulence, and abundance self- 
consistently with the MC code. We use collision rates 
calculated by Green (1992) for CS, by Flower & Launay 
(1985) for CO, by Green (1991) for H2CO, by Flower 
(1999) for HCO+, and by Turner (1995) for N2H+. Once 
the MC code calculates each level population, we simu- 
late specific molecular line profiles emitted by the model 
core at the assumed distance, d, over the spectral line 
velocities at any impact parameter by calculating the ra- 
diative transfer along rays and by convolving to match 
the spectral and spatial resolution of the observations 
with a virtual telescope simulation program. For this 
stiidy, we adopt the distance of Taurus, 140 pc, and spa- 
tial resolutions of several telescopes that we have used 
for studying star- forming cores. The simulation program 
can model molecular lines with multiple transitions and 
multiple impact parameters simultaneously. 

We model lines of CO, C^^O, CS, HCO+, and N2H+ 
at submillimeter wavelengths. Those lines have been 
observed with the 10-m telescope at the Caltech Sub- 
millimeter Observatory, the 15-m James Clerk Maxwell 
Telescope, and the 45-m telescope at the Nobeyama Ra- 
dio Observatory toward PPCs, Class and I objects (Lee 
et al. 2003; Lee et al. 2004b). We also model the ortho- 
H2CO 6 cm line, which has been observed at Arecibo 
Observatory (K. Young et al. 2004). Therefore, the sim- 
ulation of molecular line profiles and comparison with the 
actual data in different evolutionary stages will provide 
a possible opportunity to test chemistry as a tracer of 
star formation. The comparisons will appear in separate 
papers, as mentioned above. 

2.6. Characteristic Timescales 

We compare characteristic timescales for this model in 
Table 4 for a density of 10^ cm~^. We assume steady 
sate for the calculation of dust temperature, gas temper- 
ature, and the level population of molecules. The main 
heating source of the dust energetics is the radiative heat- 
ing. The energy diffusion timcscalc inside grains is much 
sniallcir than the photon absorption tiniesc;ale, i.e., the 
interval between consecutive arrival of photons. There- 
fore, the timescale of the dust energetics is mainly depen- 
dent on the photon absorption timescale, which is very 
short compared to dynamical and chemical timescales. 
In this work, we neglect small, transiently heated grains 



because they are not important to the processes that we 
consider in well-shielded regions, such as dense molecu- 
lar cores (Bernard et al. 1993). The timescale for the 
gas temperature to equilibrate with the dust tempera- 
ture through dust-gas collisions is also short enough to 
be ignored when we consider chemical and dynamical 
evolution in the high density regions. The timescale for 
molecular excitation is also much shorter than chemical 
and dynamical evolution (see Table 4), supporting the 
steady state approximation for the excitation of molecu- 
lar levels. 

We have to consider dynamics and chemistry simul- 
taneously once collapse begins {t > 0) because their 
timescales are similar within the infall radius. The exact 
process during the PPC stage is still a matter of debate, 
and cores may coalesce either through dissipation of tur- 
bulence (Myers & Lazarian 1998) or ambi-polar diffusion 
(Ciolek & Mouschovias 1993). In this work, for t < 0, 
we assume a sequence of Bonnor-Ebert spheres with a 
total timescale fixed at 10^ years with shorter timescales 
for denser stages. Although some depletion is expected 
during the early stages, the largest effect on the chem- 
istry will occur at the highest densities. Therefore, the 
relevant timescale for the PPC stage is the time spent at 
the greatest density prior to collapse. The evaporation 
timescale depends exponentially on Tp. This timescale is 
cither extremely long (To < T,.vap) and evaporation can 
be ignored, or it is extremely fast (T^ > Tevap) and can 
be considered to be instantaneous. The evaporation pro- 
cess yields molecules that drive gas phase chemistry with 
its longer timescale. From Table 4 and this discussion, it 
is clear that we need to treat the dynamical and chem- 
ical evolution together in a time dependent way, while 
we can treat the dust and gas energetics and molecular 
excitation in the steady state approximation. 

3. RESULTS OF CHEMICAL MODELS 

3.1. Standard Model 

The results of our standard chemical model, where the 
ISRF is attenuated by Ay.o = 0.5 mag, and silicate bind- 
ing energies are used for gas-grain reactions, are shown 
in Figures 6—11. In the following wc will separate the 
terms freeze-out and depletion which are commonly used 
to reference similar effects. Here we use the term freeze- 
out to rcfca' to the removal of molecules from the gas and 
onto grain surfaces. Depletion is used in reference to the 
decrease in abundance of molecular ions, which are be- 
lieved to not freeze-out directly onto grains. Ions will 
decrease in abundance, or deplete, when a parent species 
freezes onto grain surfaces (e.g. CO for HCO+). 

3.1.1. Time Dependance 

Figure 6 shows the evolution of molecular abundances 
in gas at given radii. The timescale on the x-axis starts 
from —5 X 10''', when CO reaches around its normal abim- 
dancc of 10~^. We will concentrate at first on the evo- 
lution in the central r=0.001 pc. The abundances of 
molecules increase at very early times, and decrease as 
freeze-out/depletion acts at t < 0; then at t > 0, heating 
by the central source increases the dust temperature to 
the evaporation temperatures of molecules. "'^ 

^ The evaporation temperature (Te„ap) is Tu when the evap- 
oration timescale {Tevap) is equal to the freeze-out timescale 
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After 50,000 years the dust temperature reaches the 
CO evaporation temperature (~ 25 K), so the CO abun- 
dance quickly increases. The sharp increase of CO abun- 
dance also enhances the formation rates of CS, H2CO, 
and HCN. This is more readily shown in Figure 7 where 
the time dependence of CS and H2CO is highlighted. 
Here the CS and H2CO gas abundances show the increase 
at ~50,000 years, but the ice abundance is still increas- 
ing. In the gas the formation of CS, H2CO, and HCN is 
then indirectly powered by reactions of evaporating CO 
with and He"*" . The He+ reaction in particular creates 
the fuel to power other chemistries. The CS chemistry is 
aided by the similar evaporation temperature of CO and 
sulfur atoms. 

As time increases the increasing dust temperature ap- 
proaches the evaporation temperature of the more tightly 
bound CS, H2CO, and HCN. Upon evaporation the gas- 
phase abundance of HCN and H2CO is higher than can 
be sustained via chemical equilibrium and the abun- 
dances drop. The primary process that limits the abun- 
dances of newly evaporated species is reactions with H;^ , 
and to a lesser extent He"*". Thus the presence of dust 
temperatures above the evaporation point does not neces- 
sarily dictate high abundances. However, CS shows a dif- 
ferent evolution as after evaporation its abundance does 
not drop. This is due to the freeze-out of atomic oxy- 
gen and water at small radii. If they are not frozen onto 
grain surfaces, the sulfur goes to SO over CS at late times 
in the chemical evolution. However, they exhibit signifi- 
cant freeze-out in the central region, so CS stays as the 
main sulfur-bearing species. As a result, the CS abun- 
dance at the later time becomes greater than the normal 
abundance. 

The abundances of HCO+, NH3 and N2H+ deserve 
special attention. The chemistry of HCO+ is rather sim- 
ple as it tends to follow that of CO. NH3 and N2H+ are 
also influenced by CO but instead show abundance in- 
creases as CO freezes onto grains. Thus, at r= 0.001 
pc, all molecules except for NH3 and N2H+ deplete or 
freeze-out until 50,000 years. The sharp peak of NH3 at 
this time is due to its evaporation. Similar to the other 
species the gas phase equilibrium cannot sustain the high 
abundances from the evaporated NH3 and its abundance 
drops. In the case of N2H+, CO is a major destroyer and 
the abundance of N2H+ decreases sharply in reaction to 
the CO abundance increase around 50,000 years. 

At large radii, the dust temperature reaches the evap- 
oration temperature at timescales later than at smaller 
radii, so the duration of freeze-out is longer. However, 
the amount of freeze-out is smaller because of the lower 
density, thus longer freeze-out timescale. The evapora- 
tion peaks of H2CO, HCN, and NH3 disappear at large 
radii because the dust temperature never reaches their 
evaporation temperatures. After CO evaporates, N2H+ 
is more abundant at large radii due to the higher abun- 
dance of at lower densities. 



3.1.2. Radial Dependance 

In Figure 8 and 9 we present the radial abundance 
profiles of CO and other molecules in the gas phase 

irfreeze — out ) ■ 



at selected times. ^ The left boxes of the Figures rep- 
resent abundances in gas before collapse begins. The 
CO abundance (Figure 8a) in the gas phase increases at 
t < —2.5 X 10^ years. At later times we see the signifi- 
cant effect of freeze-out which has a strong dependence 
on the radius because the freeze-out timescale increases 
going outward. When, however. To reaches the CO evap- 
oration temperature, we see a large abundance rise in the 
center. Our chemical model includes only one molecular 
species, H2 at the beginning of evolution. CO reaches the 
abundances of 10"^ and lO""^ at -9.8x 10^ and -4.7x 10^ 
years, respectively, at the center. Once Tjj reaches Tevap 
of CO (around t = 50, 000 years at r = 0.001 pc), almost 
all CO in the ice phase evaporates, so a sharp evapora- 
tion front is produced. This evaporation front propagates 
outward as the dust temperature grows with time. The 
gas-phase CO abundance increases gradually with radius 
to make a dip, where the CO freeze-out is the most sig- 
nificant, just behind the sharp evaporation front. This is 
because the freeze-out timescale increases with radius as 
shown in Figure 10. For r > 0.08 pc (Ay < 1), photodes- 
orption will remove CO from grain mantles. As the result 
of the CO chemistry depending on the dust temperature 
and density, one can simplify the abundance profile to 
three zones; evaporation, freeze-out, and normal abun- 
dance. J0rgensen et al. (2004b) used a simple "drop" 
function, which reasonably approximates this profile, to 
simulate C^^O lines to fit actual observations. 

The evolution of the abundance profile of H^, which 
initiates the chemical pathways of other molecules, is 
shown in Figure 8b. It decreases toward the center in 
all time steps because of the increase of density toward 
the center. The Hj^ abundance decreases sharply at radii 
smaller than the CO evaporation front because H;^ re- 
acts with the evaporated CO. In Figure 8c, the thick dot 
- short dash and thick solid lines compare the abundance 
profiles of CO and It^ at t= 10^ years. The long dashed 
line indicates the abundance profile if the H;^ density 
is constant through the core. This is a usual condition 
assumed in molecular clouds (Geballe et al. 1999). As 
shown, however, only at radii smaller than the CO evap- 
oration front do the model results follow the condition 
of a constant HJ density. The Hj abundance profile is 
significantly affected by the distribution of CO freeze-out 
at radii greater than the CO evaporation front. In ad- 
dition, the H;^ abundance increases with time at almost 
all radii. This results from the combination of dynamical 
and chemical evolution. Once the gas and dust are within 
the infall radius the higher abundances established in 
lower densities are carried inward. For instance, the peak 
at r= 0.08 pc stays at the same position at t= 10^ and 
2 X 10^ years because the infall radius has not reached 
the position. However, at t= 5 x 10^ years, the infall 
radius is greater than the radius of the position, so the 
peak moves inward to r~ 0.05 pc. This effect is shown 
better in Figure 9 if we compare abundance profiles at 
t= 2 x 10^ and 5 x 10^ years. 

CS (Figure 9) shows the same trend as CO, i.e., CS 
is removed from the gas phase more with time until To 
reaches the Tg^ap of CO. After To reaches Tgyap of CS, 

^ Movies of the evolution of abundance profiles as well as density, 
temperature, and molecular line profiles for all time steps are found 
in|http: / /peggysue. as. utexas.edu/sf/ CHEM/ 
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the abundance profiles are more complicated than those 
of CO. As a general trend, there arc two jumps in the CS 
profiles (shown better in Figure 10) in the inner region, 
where dust-gas chemistry is important. One is at the ra- 
dius of Tevap of CS, and the other is coincident with the 
evaporation front of CO. The first jump can be under- 
stood in the same way as CO, but the second is related to 
the CO evaporation as mentioned in the previous section. 
At radii smaller than the first peak, the abundance stays 
constant, greater than the unfrozen, normal abundance 
at outer radii. At radii greater than the second peak, the 
freeze-out dip follows, and the CS abundance increases 
to reach its normal abundance. However, at radii greater 
than around 0.02 pc, CS is dissociated by UV photons. 

H2CO (Figure 9) shows similarities to both CS and 
CO, and is frozen on grain surfaces before collapse be- 
gins. Once Td reaches its evaporation temperature, the 
abundance profile is very complex. Figure 10 shows the 
H2CO abundance profiles in two time steps, and Tgyap 
and Tfreeze-out in the previous time step. There are two 
peaks in the abundance profile. The biggest peak around 
0.001 pc is due to the evaporation of the H2CO ice built 
up in the earlier stage. At radii smaller than 0.001 pc, 
the H2CO abundance decreases toward the center, which 
is due to the equilibrium gas phase chemistry being inca- 
pable of sustaining the high abundances from the evap- 
orated ice. The peak aroimd 0.003 pc, overlaps the CO 
evaporation front, and is related to the sharp increase of 
the CO abundance. At radii greater than 0.003 pc, the 
abundance increases with radius because the freeze-out 
timescale increases with radius, so H2CO is less frozen- 
out. 

HCN shows a trend similar to that of H2CO, which has 
two peaks at the radii of CO and H2CO evaporation tem- 
peratures. However, the HCN profile shows even more 
structure. At any given time HCN shows three abun- 
dance peaks. One is due to CO evaporation and the other 
HCN evaporation. The third is related to the release of 
CN from grain mantles (CN has an estimated binding 
energy between those of HCN and CO). Gas-phase re- 
actions between Hj^ and CN result in the production of 
HCN. At 50,000 years, where there is the transition be- 
tween the first hydrostatic cores and Class objects, the 
abundance of HCN increases toward the center by about 
a factor of 10 so that the abundance at all inner radii is 
comparable to the normal abundance in spite of struc- 
ture. This is due to the evaporation of CO reaching its 
peak where all CO is returned to the gas. At this time 
reactions of CO with He+ ions reach the peak destruc- 
tive power and a small amount of carbon is transferred 
to HCN. Therefore, HCN could be a good probe to dis- 
tinguish the Class stage from the first hydrostatic core 
stage. 

NH3 is a molecule formed at later evolutionary stages, 
so it becomes more abundant with time for t < 0. How- 
ever, it also freezes out in the center before the dust 
temperature reaches its evaporation temperature, which 
is assumed to be slightly lower than that of CO. Figure 
10 shows that the evaporation front of NH3 is located 
just behind the CO evaporation front. This small dif- 
ference of the assumed evaporation temperatures of two 
molecules causes a peak between two fronts and sharp 
drop of the NH3 abundance at smaller radii than the CO 



evaporation front as the equilibrium cannot maintain the 
high abundance of recently evaporated NH3. 

HCO"*" (Figure 9) is formed by the reaction between 
CO and Ht, so it tends to follow CO. However, unlike 
CO, HCO^ shows a decrease at radii smaller than the 
CO evaporation front. This trend is related to the evo- 
lution of H^ which decreases in abundance at the small 
radii (for a given time) because of the high densities. 
In addition, HCO+ increases at all radii with time be- 
cause the infalling material carries the and HCO"*" 
abundances established in lower density outer regions to 
smaller radii. 

N2H+ (Figure 9) is formed from the interaction be- 
tween Hg" and N2, which is a molecule formed at later 
evolutionary stages and has a lower binding energy, so 
it becomes abundant later and depleted less than other 
molecules such as CO, CS, H2CO, and HCN. N2H+ is 
mainly destroyed by CO in dense regions so that it has 
an abundance hole inside the CO evaporation front at 
t > 0. This relation between CO and N2H"'" has been dis- 
cussed by J0rgcnsen ct al. (2004a) and Lee et al. (2004a). 
The abundance of N2H+ increases at all radii, except for 
r> 0.01 pc at t= 5 X 10^ years in Figure 9, with time, 
like HCO+. 

In Figure 9, the abundance of all molecules drop at 
radii greater than around 0.01 pc from 200,000 years (ma- 
genta) to 500,000 years (orange). The drop is caused by 
gas with a different chemical state being carried inward. 
The abundance peaks around 0.03 pc do not move un- 
til t= 2 X 10^ years because the infall radius is smaller 
than 0.03 pc at this time. However, the infall radius is 
^0.11 pc at t= 5 X 10'' years, so the peaks move inward 
to r~ 0.01 pc. At this moment, the CO evaporation 
front is also around 0.01 pc producing sharp drops of the 
abundances of HCN, NH3, and N2H+. 

Figure 11 compares abundance profiles of seven differ- 
ent molecules at a given time (t = 100, 000 j^ears) to il- 
lustrate the importance of the CO chemistry. There are 
two big bumps in the profiles of CS, H2CO, and HCN 
(for the moment ignoring the rise due to CN evapora- 
tion) The first one is caused by the evaporation of the 
molecules from grain surfaces. The second one is located 
at the same radius as the CO evaporation front. The 
chemistry of those molecules is affected by the sharp in- 
crease of the CO abundance. The abundance profiles of 
NH3, N2H"'", and HCO"*" are mainly dependent on that 
of CO because CO is the major destroyer of NH3 (in- 
directly) and N2H''" (directly), but the main supplier of 
HCO+. 

Rawlings et al. (1992) used the inside-out collapse 
model to examine the chemical evolution at every point 
in the envelope. They included freeze-out without 
considering desorption to conclude that molecular ions 
should be good tracers of collapsing envelopes because 
the ion destruction rate decreases as the neutral species 
are removed from the gas phase onto grain surfaces. 
However, Choi et al. (1995) did not find any evidence 
for the freeze-out of CS in B335. Our models with des- 
orption can explain the abundant CS during protostellar 
collapse {t > 0) as we described above and in §3.1.1. 
In addition, based on the narrowness of emission lines 
in B335 (Mcnten ct al. 1984), Rawlings et al. (1992) 
suggested that the NH3 lines could not trace high veloc- 
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ity infalling material due to the freeze-out. However, we 
have found that the drop of the NH3 abundance, at small 
radii at f > 0, is caused from the CO evaporation rather 
than the freeze-out of NH3 itself. Therefore, desorption 
should not be ignored even in low-mass star-forming re- 
gions. 

We note that we have not included any grain surface 
chemistry beyond the production of gas-phase H2 and 
H2O ice. Some molecules such as NH3, H2CO, CH4, and 
CH3OH could potentially form via grain surface reactions 
and, if so, their abundances in the inner "hot core" could 
be enhanced in the immediate evaporation states. 

3.2. Other Models 

Figure 12 and 13 show abundance profiles modeled 
with three different external visual extinctions, at the 
initiation of collapse {t = 0) and t = 100, 000 years. Note 
that the edge of the cloud in Figure 12 and 13 is not di- 
rectly exposed to the ISRF due to Av,e- The main differ- 
ences are higher abundances at the outer radii for higher 
Av,e, caused by lower photodissociation of molecules by 
the interstellar UV photons or lower dissociative recom- 
bination with electrons. However, CO does not show any 
difference among different Ay.c because CO photodisso- 
ciation is not included. CS is less abundant even at inner 
radii as well as outer radii in the standard model with 
Av,o = 0.5 mag than in the models of higher Av,c at 
t = (Figure 12). This is due to the low extinction in the 
early, pre- freeze-out, Bonnor-Ebert evolutionary phases, 
which favors CS formation over SO. At such early times 
models with higher extinction produce more SO which 
freezes onto grain surfaces. In the standard model, CS 
is more abundant at the radii smaller than 0.0025 pc 
at 100,000 years (Figure 13). The abundance difference 
at inner radii is the biggest between 0.001 pc and 0.002 
pc. In the standard model, the CS abundance does not 
change much (by a factor of 2), but it varies by the 2 
orders of magnitude in the model with Av,o — 3 mag. 
This large effect is due to the difference in the overall 
evolution of sulfur pool. In the standard model there 
is higher amount of S atoms frozen on grains which is 
seen in a reduction of gas-phase CS abundance in Figure 
12. Sulphur atoms are assumed to evaporate at similar 
temperatures to CO and thus the sudden return of sul- 
fur and carbon to the gas fuels CS production. For the 
higher extinction cases there is a large amount of SO on 
grains which when evaporated does not lead to efficient 
CS creation. The CS abundance at the cloud edge for the 
model with Av,c = 0.5 mag does not show a rapid de- 
cline from photodissociation due to a higher abundance 
of ionized carbon at r=0.01 pc which results in efficient 
CS formation. The ionized carbon exists due to slower 
chemical timescales in the outer layers. 

Figure 14 and 15 compare results of models with differ- 
ent binding energies in dust grain reactions. The bind- 
ing energies to the CO mantle and the H2O mantle are 
smaller and greater than that to the Si02 mantle by the 
ratio of 0.82 and 1.47, respectively. At t = 0, the model 
with higher binding energy exhibits greater freeze-out of 
CO, CS, and HCN onto dust grains and greater deple- 
tion of HCO+. Even NH3 and N2H+ show significant 
freeze-out /depletion in the model with higher binding 
energy at this early evolutionary stage. The results of 
the model with the higher binding energy is much closer 



to the standard model results for CO, CS, HCN, and 
HCO+, while the opposite is seen for NH3 and N2H+. 
Another difference is that HCN shows little freeze-out in 
the model with low binding energy. Therefore, we can 
rule out some dust properties from the PPC observations 
by comparing HCN with NH3 and N2H"'". For example, if 
no freeze-out/dcpletion is seen in these three molecules, 
dust grains are covered by CO. On the other hand, dust 
grains have water dominant mantles if all three molecules 
show freeze-out /depletion. If only HCN shows freeze-out, 
bare Si02 grains are the main dust component. 

After collapse begins (t = 10'"' yr case), the dust tem- 
perature increases, and the difference between various 
binding energies is more significant and complicated. 
Greater freeze-out of CO, CS, and HCO+ depletion is 
observed in the model with the higher binding energy. In 
addition, the evaporation fronts of molecules are located 
at smaller radii due to higher evaporation temperatures. 
In contrast to the case at t=0, the results of the model 
with lower binding energy are closer to those of the stan- 
dard model in CO, CS, HCN, and HCO+. Sharp differ- 
ences are also seen for NH3 and N2H+ where the model 
with high binding energy shows significant abundance 
structure not seen in the other cases. This structure is 
caused by the freeze-out of N2 onto H2O mantles. Thus 
the abundance structure for NH3 and N2H+ is due to the 
evaporation of N2 from the grains. Models with lower 
binding energy have little N2 freeze-out and thus both 
NH3 and N2H"'" are not directly affected by N2 evapora- 
tion. 

Bergin & Langer (1997) used the same chemical code 
to study the chemistry in the PPC stage. They showed 
that CO and HCO+ were not depleted on CO grain sur- 
faces, unlike our result. This discrepancy is caused by the 
different density evolution. In their model, the density 
does not increase above 10^ cm~^, but it goes up to 10^ 
cm~^ at the center in our model. The higher density re- 
duces the freeze-out timescales of CO, so CO freezes-out 
even on CO dominant grain mantles, and HCO+ depletes 
from the gas phase. 

4. RESULTS OF LINE SIMULATIONS 

We use a 30" beam, which is the typical beam size of 
the CSO 10 m telescope around 230 GHz, for CO 2 - 1, 
and HCO+ 3-2 lines. The CS 2 - 1 and HCO+ 1 - 
lines have also been simulated with a 30" beam. For the 
N2H+ 1 - fine, the beam size (17") of the NRO 45 m 
telescope at 93 GHz has been used. The beam sizes for 
the C^^O lines have been adopted from J0rgensen et al. 
(2002); 33", 21", and 14" for 1 - 0, 2 - 1, and 3-2, 
respectively. We assume that the ortho-H2CO 6 cm line 
is observed with the Arecibo telescope {di, = 60"). The 
HCN 1 — line is also simulated with various resolutions. 
We use the same velocity resolution of 0.1 km s~^ and 
the same microturbulence of 0.15 km for all lines. 

Figure 16 shows the evolution of the CO 2 — 1 and 
the CS 2—1 lines. In our models the exclusion of CO 
photodissociation may modify the results with respect to 
line profiles. However, the additional extinction included 
in the model limits these effects. The CO 2—1 line of the 
model with Av,e = 3 mag is much weaker than that of 
the model with Ay.o = 0.5 mag. The CO 2—1 transition 
has low critical density, so it can be thermalized even at 
the outer radii, and the kinetic temperature in the for- 
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mer is much lower than that in the latter at large radii 
as shown in Figure 5. In addition, The CO 2 — 1 line 
of the model with Av,e = 3 mag shows the blue asym- 
metry caused by infall. However, The CO 2—1 line of 
the model with Ay.e = 0.5 mag does not change much in 
all time steps except in the time step of 500,000 years. 
The line at 500,000 years shows the inverted asymme- 
try rather than the blue asymmetry. This is caused by 
the distribution of its excitation temperature. The ex- 
citation temperature in the model with Ay.e = 0.5 mag 
increases at radii greater than 0.016 pc (see Figure 17), 
causing the inverted asymmetry. This result suggests 
that we have to consider simultaneously parameters such 
as proper kinetic temperature and abundance profiles as 
well as density and velocity structures. 

The CS 2 — 1 line does not show a big difference in 
intensity between two models while the blue asymmetry 
profile is more prominent in the model with Av,e = 3 
mag. The similarity in intensity may be caused by the 
higher critical density of about 10^ cm~'^ of the line, so 
lines are emitted dominantly from the inner, denser re- 
gion, where the excitation temperatures of two models 
are the same. However, the more prominent blue asym- 
metry in the model with Av,e = 3 mag can be caused 
by the higher CS abundance in the outer part in the 
model with the higher Av,e- The wing of the line be- 
comes broader with time because the infall velocity in- 
creases with time. Since the material of the model core 
falls onto the central object, the density within the in- 
fall radius decreases with time, and the peak radiation 
temperature of the CS 2—1 line becomes smaller. 

Figure 18 illustrates the importance of the CO chem- 
istry, which causes a peak in the CS abundance profile, 
even in line profiles. If the second bump of the CS abun- 
dance profile at 100,000 years, which is caused by the 
CO evaporation, is ignored, then the CS 2 — 1 becomes 
weaker, and the wing structure disappears. Therefore, 
the wing structure in the CS 2 — 1 line profile could be 
evidence that the dust temperature is high enough to 
evaporate CO or CS. However, this wing structure could 
be confused by an outflow, which is not considered in 
this work, so high spatial resolution maps are important. 

Figure 19 shows the evolution of HCO+ 1 — and 3 — 2 
lines. The lines become broader with time due to greater 
infall velocity at later times. The blue asymmetry pro- 
file, which is an indicator of the infall motion, is robust 
in both models with two different external extinction en- 
vironments. In addition, the blue asymmetry of the 1 — 
and 3 — 2 lines grow with time, as seen in Gregersen ct 
al. (1997). Both 1 - and 3 - 2 lines of the standard 
model are weaker than those of the model with Av,e = 3 
mag at 500,000 years. This is caused by the dissociative 
recombination of HCO+ with electrons at large radii in 
the model with the lower Ay.c- 

The evolution of the ortho-H2CO 6 cm line is shown in 
Figure 20. The absorption in that line occurs against 
the cosmic background radiation (CBR) because of a 
coUisional pumping that lowers the excitation temper- 
ature below TcBR (Towncs and Cheung 1969, Zhou et 
al. 1990). The pumping works up to about 10^ cm~^ in 
density, so the line goes into emission at sufficiently high 
densities if the abundance of H2CO is constant through 
a entire core. However, the abundance profiles of H2CO 
are complicated with many variations. As shown in Fig- 



ure 9, H2CO shows significant freeze-out before the dust 
temperature reaches the CO evaporation temperature. 
In addition to the variation of H2CO with time, its abun- 
dance profile varies with Av,e, that is, the cloud environ- 
ment. Therefore, the evolution of the line depends on 
the external visual excitation. In a lower Ay.c environ- 
ment, the abundance drops sharply at outer radii, where 
the line cannot be thermalized (Figure 26). In that re- 
gion, the pumping mechanism works efficiently due to the 
low abundance and the high kinetic temperature so that 
deeper absorptions are produced. On the other hand, at 
higher Av,c, H2CO is protected against the photodisso- 
ciation, and the lower kinetic temperature reduces the 
collisional pumping to produce an emission line at 6 cm. 
Based on this result, we can qualitatively distinguish the 
different environments of a star-forming core; an emis- 
sion line at higher Ay.c and an absorption line at lower 
Ay^e- We see only absorption lines toward the centers of 
three PPCs in K. Young et al. (2004), which suggests 
that the cores are embedded in relatively low Ay^g- This 
is consistent with the higher Go that they found from the 
CO observations, compared to what is expected from the 
Ay^e found in the dust radiative transfer calculation. 

Figure 21 shows the evolution of the N2H+ 1 — iso- 
lated hyperfine component (Iqi — O12). The line has a hy- 
perfine structure with 7-components, with 6 components 
located close to each other. The N2H+ abundance in- 
creases until CO evaporates from grain surfaces. There- 
fore, N2H+ can be a good tracer of the density structure 
in the earlier evolutionary stages before To reaches T^yap 
of CO. Once CO starts to evaporate {t = 50, 000 years), 
the N2H+ line becomes weaker. N2H+ is abundant even 
at radii larger than 0.03 pc in the model with Ay^e = 3 
mag (see Figure 12), so the line becomes optically thick 
and shows the self-absorption profiles even in earlier time 
steps. However, the blue asymmetry in this line traces 
the infall at the radii outside the CO evaporation front 
because the N2H+ abundance drops by a few orders of 
magnitude inside the CO evaporation front. 

5. DISCUSSION 

In actual observations, we have to choose a proper 
beam size, transition, and observational technique such 
as central position observation or mapping in order to get 
correct information about star formation. HCN could be 
a good probe for distinguishing Class from the FHSC, 
as described above. However, we have to use a very good 
resolution for the classification. Figure 22 shows the dif- 
ference between simulated observations of the HCN 1-0 
line with two beam sizes of 5" and 50". The HCN 1-0 
has a hyperfine structure of 3-components, which are well 
separated from each other. In the observation with a 
lower resolution, it is difficult to tell when the transi- 
tion between the FHSC and Class occurs. On the 
other hand, a big change in the line strength between 
two stages is seen in the simulation with a high resolu- 
tion. Therefore, an interferometer observation is needed 
to differentiate Class from the FHSC with HCN fines. 

The chemical evolution is affected by the environment 
of cores as shown in Figure 12 and 13, as well as grain 
properties inside cores as shown in Figure 14 and 15. The 
environment of star-forming cores can be described by 
the combination of the external visual extinction and the 
strength of the UV field. The UV field is extincted dif- 
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ferently by different grain properties along lines of sight, 
such as the material and size of dust grains, and by differ- 
ent cloud geometries. Av,e is important in the chemistry, 
and the strength of the UV field is important in the gas 
energetics. 

In the following we compare how various environments 
affect molecular line profiles. Figure 23 shows the com- 
parison of the C^*0 1 — 0, 2—1, and 3—2 lines in 
the models of two different visual external extinctions 
(Av,e = 0.5 and 3 mag) at 200,000 years. Here, Go is 
calculated with the equation shown in §2.3 assuming the 
same grain properties and cloud geometry. All linos are 
simiilatod toward the center of the model core. Wo use 
the same beam size for each line as those of j0rgensen et 
al. (2002); 33", 21", and 14" for 1 - 0, 2 - 1, and 3-2, 
respectively. The 1 — line has similar strength to the 
2—1 lino, which agrees with the actual observations and 
the simulated linos with an abundance profile of a drop 
function (J0rgensen et al. 2004b). The wing structure 
in the C^^O lines is due to the infall velocity structure 
rather than an optical depth effect. The optical depths 
of the lines are <C 1. The 1 — line does not change with 
different Av,c, but the 2 — 1 and 3 — 2 lines are weaker 
for higher Av,e- Therefore, multi-transition observations 
can be good probes of the cloud environment. 

Figure 24 compares the ortho-H2CO 6 cm lines along 
various lines of sight at 50,000 years to show the effect 
of different Av.c- For smaller Av,c, H2CO has a sharp 
drop in the abundance profile duo to the photodissocia- 
tion, and Tk is high because of the PE heating at outer 
radii as shown Figure 26. Therefore, the excitation tem- 
perature of the 6 cm line for the model with the smaller 
Av,e is much lower than the temperature of the CBR 
{TcBR = 2.7 K) at larger radii than 0.025 pc. The 6 cm 
line produces deep absorption at small impact parame- 
ters, but the absorption dip is much weaker at a large 
impact parameter because of the low abundance at large 
radii. On the other hand, in the core with greater Av,e) 
even though the H2CO abundance is very high at outer 
radii, the absorption dip is not deeper than at inner radii 
because the excitation temperature is similar to 2.7 K. 
As a result, the 6 cm line goes into emission at the center. 

We show the effect of different UV strengths in Figure 
25. We use the same Ay.c of 3 mag, but two different Go 
have been considered. One is obtained from the equa- 
tion in §2.3 (Go = 0.005), and the other is assumed to be 
10 times bigger (Go = 0.05). The abundance profiles of 
two models are almost the same because the chemistry 
is not affected much by the kinetic temperature at later 
time steps as shown in Figure 26. However, the kinetic 
temperature of the model with Go = 0.05 is higher than 
that of the model with Go = 0.005 at outer radii. In 
the model of Go = 0.05, the excitation temperature of 
the 6 cm line is below Tcbr at larger radii than 0.04 
pc. Therefore, the high abundance and the lower ex- 
citation temperature than 2.7 K at outer radii produce 
deep absorption lines at all impact parameters. These 
comparisons suggest that we need to do multi-position 
observations as well as multi-transition observations in 
order to test the cloud environment. 

We have to note that there is a caveat about collision 
rates. Collision rates, in theory, have been calculated 
down to 10 K, but, in our model, Tk goes down to about 
6 K at the center in the PPC stage and the FHSC stage. 



Wc, therefore, extrapolate the collision rates linearly to 5 
K. Wc also emphasize that outflows and hot cores, where 
the dust temperature is higher than 100 K, are not con- 
sidered in this work, and they can affect the interpreta- 
tion of lines from the protostellar stage. 

6. CONCLUSIONS 

Although we still have limitations such as 1-D geom- 
etry and no grain surface chemistry in this proposed 
model, we certainly see that the chemical timescale is not 
short or long enough to ignore the abundance variations 
of species during the dynamical evolution in star forma- 
tion, especially in the inner dense regions as shown in Ta- 
ble 4. The abundances of molecules vary by a few orders 
of magnitude with time at the central part, but the vari- 
ation becomes less significant in outer regions. There- 
fore, abundance profiles at given times arc very wiggly. 
Even though there are small bumps and wiggles in the 
abundance profiles of molecules, the overall structure of 
the profiles can be understood as interactions between 
freeze-out and evaporation of the molecules, themselves, 
and CO. Some of the bumps, jumps, or drops are due to 
CO evaporation. The presence or absence of CO from the 
gas phase has a tremendous effect on the less abundant 
species. 

Astronomers use molecular lines to study star-forming 
regions. However, the lino observations can be misun- 
derstood if one docs not consider the abundance varia- 
tions with time or with radius, which have been shown in 
this study. According to our results, carbon- or sulfur- 
bearing molecules such as CO, CS, and H2CO readily 
freeze-out from the gas on all kinds of dust grains as 
shown in other theoretical modeling work. HCO^, a 
daughter species of CO, also disappears from the gas 
phase due to the CO freeze-out. However, nitrogen- 
bearing molecules such as NH3 and N2H+ could be a 
good probe of the density structure in the PPC stage if 
grain mantles are covered with Si02 or CO dominantly. 
On the other hand, they are frozen-out/depleted signif- 
icantly on water dominant grain mantles at very late 
times of the PPC stage or at the FHSC stage. HON 
is an interesting species. It exhibits strong freeze-out 
on silicate or water dominant mantles as other carbon- 
bearing species, but it becomes abundant with time on 
CO dominant grain mantles at t < 0. Therefore, we may 
use HCN with NH3 or N2H+ to learn about grain proper- 
ties in PPCs. In addition, HCN becomes abundant very 
quickly in the transition from the FHSC to Class 0, so 
it could be a probe to distinguish the two evolutionary 
stages. In the later stages than the FHSC, the assump- 
tion of a constant abundance or a simple functional form 
of the abundance structure in some molecules can lead 
one to incorrect density and velocity structures because 
abundance profiles of those molecules vary significantly 
with time and radius, so they are very wiggly, especially, 
in the very inner region. We, therefore, suggest that the 
chemical evolution should be considered simultaneously 
with the dynamical evolution when one studies the evo- 
lution in star formation, and high spatial resolution maps 
are very important. However, one can still adopt simpli- 
fied functional forms of abundance profiles depending on 
the types of molecules. Molecules such as CO, HCO"*", 
NH3, and N2H+ show rather simple abundance struc- 
tures that can be approximated by a step function or a 
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drop function. 

In our study wc find strong evidence for line wings pro- 
duced in the infalling envelopes. The primary reason for 
the presence of the infall wings is the sharp rise of molec- 
ular abundances within the specific evaporation front for 
that species. However, these same species can likely be 
liberated from grains within the outflow, which can hide 
the presence of wings produced solely from infall. 

The complicated abundance profiles of HCN, CS, and 
H2CO highlight another important result of this work. 
To first-order, abundance profiles of species in star- 
forming cores can be approximated using simple three 
component profiles (high in center, freeze-out in middle, 
and an abundance rise in outer layers). However, there 
are important second-order effects where abundances in 
the gas are significantly altered by the evaporation of 
separate species. These second-order effects could be 
important if, for example, they resonate with particu- 
lar changes in the velocity field and may eventually ac- 
count for some difficulties in the interpretation of obser- 
vations. The infiuence of CO on the chemistry of CS 
is one example to show how the second-order effect is 
coupled with the velocity field in collapse as shown in 
§4. J0rgensen et al. (2004b) calculated correlations be- 
tween various molecular abundances and envelope masses 
of star-forming cores. According to their analysis, CO 
and HCO"*" show very high correlation with the envelope 
mass (i.e. evolutionary stage), but SO and HCN show 
little correlation with the envelope mass. These results 
could be explained by the second-order effects seen in 
our models. CO does not have the second-order effect, 
and the abundance of HCO"*" is directly linked to CO. 
Therefore, CO and HCO+ both become more abundant 
because of the CO evaporation as the envelope mass de- 
creases. However, HCN is affected by CO and CN evap- 
oration, and the second-order effects may complicate the 
relation between the HCN abundance and the envelope 
mass. For the SO chemistry, oxygen plays a role. Oxy- 
gen is locked onto grain surfaces in the form of water, 
and sulfur-bearing species, including SO, react with car- 
bon, which is produced by the evaporated CO and He"*" 
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ions, to form CS. As a result, the evaporation of SO does 
not result in the high abundance of SO in the gas phase, 
so a tight correlation between SO and the envelope mass 
is not seen. 

Different timescales during the PPC stage and differ- 
ent dynamics during collapse need to be explored in fu- 
ture works. In this work, a few simple tests of different 
timescales, that is, different density evolutions, during 
the PPC stage have been done. According to the re- 
sults, the freeze-out/depletion dip is shallower and nar- 
rower during the PPC stage if the timescale of the PPC 
stage is shorter than that of standard model, and vice 
versa. However, during collapse, the depth of freeze- 
out /depletion dip is the same regardless of the timescale 
of the PPC stage, and the difference between the vari- 
ous timescales becomes small with time. CS and HCN 
are the most sensitive molecules to the density evolution 
during the PPC stage. Their abundances vary more than 
one order of magnitude depending on the timescale of the 
PPC stage. 

As separate work, we are also comparing this model 
with actual observations. We have compared the result 
of H2CO with the ortho-H2CO 6 cm line toward three 
different PPCs (Young el al. 2004). Evans et al. (2004) 
studies the chemical status of B335 by comparing multi- 
transition observations of various molecules with simu- 
lated ones with chemical models of different dust prop- 
erties and initial conditions. Comparisons of the results 
of this work with Multi-transition and multi-position ob- 
servations done with single dishes and interferometers in 
more Class and I sources will appear in separate papers 
(Lee et al. 2004a, 2004b). 
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TABLE 1 

The summary of the standard model 





Ro'^ 




M "1 


Av.e" 


Go' 


pc 


pc 


km s" 1 


M0 


mag 




0.00033 


0.15 


0.22 


3.6 


0.5 


0.4 



^Inner radius of tlic core 
''Outer radius of tlic core 

Sound speed for calculating density structures in inside-out collapse 
''The initial mass of core 
"External Visual extinction 

^The strength of the UV field relative to the local ISRF 



TABLE 2 

The surface binding energy of species onto Si02 dust grain mantle 



species 



H2O 
CS 

H2CO 
HON 
CN 
CO 

NH3 



binding energy (K) 



5700 
1999 
1722 
1722 
1476 
1181 
1082 



TABLE 3 
Initial elemental abundances 





Abundance 


Element 


(Relative to H2) 


He 


0.28 


Hc+ 


1.33 X 10""' 


C+ 


1.46 X 10-" 


N 


4.28 X 10-5 





3.52 X 10-" 


Si+ 


4.0 X 10"" 


Mg+ 


6.0 X 10"" 


S+ 


4.0 X 10-* 


Fe+ 


6.0 X 10-" 


Na+ 


4.0 X 10-11 


P+ 


6.0 X 10-8 
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TABLE 4 
Characteristic Timescales 





Process 


Timescale (yrs) 


Timescale (yrs) 
for n(H2) = 105 cm-3 


Dynamics 


Bonnor-Ebert spheres (t^O) ^ 
Free-fall (t>0) 


dependent on central densities 
2.4x10'^ 


D X lU 

7.6 X 10* 


Dust energetics 


Radiative heating 


r j c ^ 


U.4<5 


Gas energetics 


Gas-dust collision 


1.26x10° 

^/rr"(H2)2 


0.004 


Chemistry 


Gas-phase chemistry 
Freeze-out 
Evaporation (desorption) 


2.6x10'' 
n(H2) 
6X10^ / m f 
"(H2) V 


2.6 X 10* 
10^ 

6.5 X 10^1 (Td = 10 K) h 
0.004 (Td = 30 K) 


Excitation of molecular 
energy levels 


Collisions 


3.2x10^ 
~ V^"(H2) 


0.013 ' 


^ For t < 0, wc have assumed 


a sequence of Boniioi-Ebert spheres 


with a total time fixed at 10® years. 


As the central density of the 



Bonnor-Ebert sphere increases, the timescale decreases. 
''We adopt the timescale for the Bonnor-Ebert sphere of ric = 10^' cm~^. 
'^r (= 0.15 pc) is the size of our model core, and c is the speed of light. 
'^Tfc = 10 K is used for all timescales. 
"Induced by the cosmic-ray ionization 

is the molecular weight. 
^Ef, is the binding energy of a molecule onto the surface of dust grains. 
^For CO 

'De-excitation timescale of the transition of CO J= 1 — 0. The de-excitation rate is adopted from Flower and Launay (1985). 
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The Sequence of Modeling 



Physical Model 
in a given Dynamics 



gas to dust 



Continuum Radiative 
Transfer 




Gas Energetics 
(heating & cooling 
>alance 






Simulation of Line 
Profiles 
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Fig. 1.— 



The flowchart to show the sequence of the whole modeUng process. 
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Fig. 2. — Density and velocity profiles of the model core with time. In the left box, dashed lines represent density profiles of Bonnor-Ebert 
spheres of nc = 10"', 3 X 10*, lO'', 3 X lO"', lO'', 3 X 10^', and 10^. The dotted hne indicates the very initial constant density profile. Solid 
lines represent density profiles during collapse at 0, 2.5 x 10^, 5 x 10^, lO'*, 2.5 x lO'*, 5 x lO**, 10^, 1.6 x 10^, and 5 x 10^; the timescale 
during collapse is divided into 190 time steps. The right box shows velocity profiles of collapse at the same time steps as in density profiles. 
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Fig. 3. — Physical parameters of given parcels of gas with time; the left upper box for positions, the right upper box for 
densities, the left lower box for dust temperatures, and the right lower box for total visual extinctions, respectively. Each line 
represents the parcel of the grid number of 50, 100, 150, 200, 250, 300, 350, 400, 450, and 500 among 512 grid points. A time of 
zero indicates the starting point of the gravitational collapse of the model core. In the position plot, a parcel stays in a given 
position until the infall wave reaches its initial position, and then moves inward with time. There is a discontinuity in density 
and Av at t = because we combine the Bonnor-Ebert spheres and inside-out collapse artificially. In the dust temperature plot, 
before the beginning of collapse, the temperature of a given parcel decreases with time because density and visual extinction 
increase, so the ISRF, the only heating source of PPCs, cannot penetrate well. However, after collapse starts, the core is heated 
by accretion, and temperatures go up quickly to produce a big jump around 50,000 years. There is another sharp increase of 
temperature caused from the turn-on of photospheric luminosity from gravitational contraction and deuterium burning in the 
central star at 100,000 years (see §2.2). X-axis of the dust temperature plot is not the same as that of other plots. The variation 
of dust temperature is small in early stages, so we plot it from t=-300,000 years to show better resolution in later stages. 
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Fig. 4. — The results of dust contimuiin radiative transfer ealculation. The left box shows dust temperature profiles of Bonnor-Ebert 
spheres of seven different central densities same as in Figure 2. The dotted line in this box represents the temperature profile for the initial 
density profile. In the right box, lines represent dust temperature profiles during collapse at the same time steps as in Figure 2. There is a 
big jump between 20,000 and 50,000 years. The model core is in the stage of the first hydrostatic core until 20,000 years. At that stage, 
accretion luminosity is still small because of a large radius (5 AU) of the central source; consequently the dust temperature stays very 
low. The dotted line and dashed line show another jump in temperature caused from the turn-on of pliotospheric luminosity at 100,000 
years. The small box inside the right box shows the evolution of the total luminosity, which is the sum of the accretion luminosity, the disk 
luminosity, and the protostellaj: luminosity, in the unit of the solar luminosity. 
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Fig. 5.- Kinetic temperatures calculated with the energetics code. The selected time steps are the same as in Figure 2 and 4. Left 
two boxes show gas temperatures in seven Boimor-Ebert spheres. The dotted lines indicate the gas temperatures calculated in the initial 
state. Gas temperatures during collapse are represented in right boxes. The upper boxes and the lower boxes compare results of heating 
by the ISRF atteimated by different magnitudes of external visual extinction. As the central density of the B.E. sphere grows, the gas 
temperature at small radii gets close to the dust temperature. 
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Fig. 6. — The evolution of abundances at given radii, r= 0.001 pc, 0.005 pc, and 0.01 pc. The lines arc for CO (solid), CS (dot), H2CO 
(short dash), HCN (long dash), NH3 (dot - short dash), HCO+ (dot -long dash), and N2H+ (short dash - long dash) from the top to the 
bottom. Lines are shifted up and down to compare with CO evolution. CS is shifted up by the 1.7 orders of magnitude, and HCN, NH3, 
HCO+, and N2H+ are shifted down by the 5, 10, 15, 18, and 21 orders of magnitude, respectively. 
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Fig. 7. — The evolution of the abundance of CS and H2CO at a given radius, r= 0.001 pc. SoUd and dotted lines represent abundances 
in gas and in ice versus time. The y-axis is not a logarithmic scale, unlike in Figure 10. 
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Fig. 8. (a) and (b). CO and abundance profiles with time in some selected time steps; the left box for the evolution of Bonnor-Ebort 
spheres and the right box for the evolution during collapse in each window. In the left box, lines arc for time steps of 5 X 10''' (thin solid), 
2.5 X 10^ (dot), 1.25 X 10^ (short dash), 6.5 x 10" (long dash), 3.5 X 10"' (dot - short dash), 1.5 X lO"* (dot - long dash), 5 X 10'^ (short dash 
- long dash) years before collapse, and the starting point of collapse (thick solid). The right box represents the CO abundance profiles at 
10^ (sohd), 10" (dot), 2 X lO" (short dash), 5 X lO" (thin long dash), 5.2 X lO" (thick long dash), 10"' (thin dot - short dash), 1.5 X 10^ 
(thick dot - short dash), 2.0 X 10''' (dot - long dash), and 5 X 10^' (short dash - long dash) years after collapse begins, (c). The effect of 
CO and dynamical evolution to the abundance profiles. The long dashed line indicates the abundance profile if the density is 
constant. The thin dotted lines show how the abundance profile established in an earlier time step affects an abundance profile in a later 
time step. 
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log r (pc) log r (pc) 

Fig. 9. — The evolution of abundance profiles. The left box for each molecule shows the molecular abundance profiles before collapse, 
and the right boxes represent the abundance profiles after collapse begins. The selected time steps in left boxes are 5 X 10^ (thin black), 
2.5 X 10^ (red), 1.25 X lO^^ (green), 6.5 X 10"* (blue), 3.5 X lO* (cyan), 1.5 X 10"* (magenta), 5 X 10^ (orange) years before collapse, and the 
thick black line indicates the abundance profile at the initiation of collapse (t=0 year). Time steps in right box are 10^ (black), 10"' (red), 
2 X lO"* (green), 5 X 10'' (blue), 10^ (cyan), 2.0 X lO'' (magenta), and 5 X 10^ (orange) years after collapse. Moves for all time steps are 
found in http://peggysue.as.utexas.edu/sf/CHEM/ 
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Fig. 10. — The comparison of evaporation tiincscalc and frcczc-out tinicscalc. Thick soUd and dashed lines represent abundance profiles 
in gas and ice at t = 100, 000 years. Thin solid and dotted lines indicate the freeze-out timescale and evaporation timescale, respectively, 
a.t t = 100, 000 years. At the point when two timescale are same, molecules are desorbed out from dust grain surfaces quickly producing a 
sharp evaporation front. The evaporation front moves outward with time. 
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Fig. 11. — Abundance profiles at a given time, t = 100,000 years. The distributions of other molecules are mainly dependent on that 
of CO. The abundances decrease at the outer radii because molecules are easily photodissociated for low external extinction (0.5 mag). 
Profiles have been shifted up and down to compare with the CO abundance profile better, so y-axis does not give correct abundances of 
molecules except for CO and HCN. CS and H2CO are shifted up by the 3.5 and 1.5 orders of magnitude, and NH3, N2H+, and HCO+ axe 
shifted down by the 2.5, 2, and 5.5 orders of magnitude, respectively. 
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Fig. 12. — Abundance profiles in models with different external visual extinctions, at the initiation of collapse. 
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Fig. 13. — The same as in Figure 12, but at t = 100, 000 years. 
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Fig. 14. — Comparison between different binding energies onto dust grain surfaces at the initiation of eoUapse. Solid line indicates the 
result of our fiducial model with the binding energy onto Si02 grain surfaces. Dotted line represents the result of a model with higher 
binding energy onto water dominant grain surfaces. Dashed line is for lower binding energy onto CO dominant grain surfaces. 
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Fig. 15. — The same as in Figure 14, but at t = 



100, 000 years. 
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Fig. 16. — The evolution of CO 2-1 and CS 2-1 lines. Solid and dotted lines represent line profiles simulated from the chemical model 
with 0.5 mag and 3 mag of external visual extinctions, respectively. Numbers in upper part of boxes indicate time steps. Line profiles at 
the earlier time steps than 10,000 years are very similar to line profiles at 10,000 years. 
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Fig. 17. — The profiles of excitation temperature of the CO 2—1 fine in models with different external visual extinctions at 500,000 
years. 
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Fig. 18. — The comparison of the CS J= 2 — 1 hnes with different abundance profiles at t=100,000 years. Sohd hne has been simulated 
with the actual abundance profile obtained from chemical model, and dotted line has been simulated with the abundance profile where 
the second bump has been removed artificially. The solid line and dotted in the small box indicate the modeled abundance profile and the 
artificial abundance profile, respectively. 
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Fig. 19. — The evolution of HCO"*" 1—0 and 3—2 lines. Solid and dotted lines represent line profiles simulated from the chemical models 
with 0.5 mag and 3 mag of external visual extinctions, respectively. 
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Fig. 20. — The evolution of H2CO 6 cm line. Solid and dotted lines represent line profiles simulated with the results of chemical models 
of Av,e = 0.5 and 3.0 mag, respectively. 
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Fig. 21.— 



The evolution of N2H+ isolated component (Iqi — 0i2). 
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Fig. 22. — 



Comparison of HCN 1 — lines simulated with different resolutions. 



38 




39 




-2-1012-2-1012-2-1012 

Vlsr (km s"') 



Fig. 24. — Comparison of ortho-H2CO 6 cm line in Av,e = 0.5 and 3 mag at 50,000 years in the assumption of the same dust properties. 
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Fig. 25. — Comparison of ortho-H2CO 6 cm line in Go = 0.05 and 0.005 with the same external visual extinction of 3 mag at 50,000 
years. 
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Fig. 26. — Models of various Av,e and Go at t = 50,000 years. In the right bottom box, lines represent the excitation temperatures of 
the strongest component in the hyperfine structure of the ortho-H2CO 6 cm line. The thick horizontal line indicates the temperature of 
the CBR (2.7 K). 



